library(ggplot2)
library(ggthemes)
library(lattice)
library(plot3D)
library(rriskDistributions)
library(stringr)
require(MASS)
## Loading required package: MASS
time.ini = Sys.time()
Fonction de lecture des stats résumées. La fonction est prévue pour lire les fichiers de stats résumées dans des répertoires de la forme “Ne100/simu_1” (en modifiant les valeurs de Ne et de simulations). Elle renvoie les résultats sous forme de data frame en ajoutant aux statistiques résumées les colonnes correspondant à l’effectif efficace et la simulation.
#Data frame construction
data.frame.stat = function(seq_ne, max_gen, max_simu){
#Créer matrice vide qui contiendra les stats
mat_stat = matrix(data = NA, nrow = length(seq_ne)*max_gen*max_simu, ncol = 64)
mat_stat = as.data.frame(mat_stat)
#Compteur pour indiquer les lignes à remplir
cpt = 0
for (ne in seq_ne) {
if (dir.exists(str_c("Ne", toString(ne), "/"))) {
dir_ne = str_c("Ne", toString(ne), "/") #Répertoire de taille effective
row_min = cpt*max_gen*max_simu+1 #Ligne min. où on écrit les stats pour le Ne donné
row_max = (cpt+1)*max_gen*max_simu #Ligne max. où on écrit les stats
mat_stat[row_min:row_max,1] = toString(ne) #Ecriture du Ne correspondant pour les stats qui seront écrites
cpt = cpt+1 #Incrémentation du compteur
for (nb_simu in 1:max_simu) {
dir_simu = str_c("simu_", toString(nb_simu), "/")
row_min_simu = (nb_simu-1)*max_gen + row_min
row_max_simu = nb_simu*max_gen + row_min -1
mat_stat[row_min_simu:row_max_simu,2] = nb_simu
file_path = str_c(dir_ne, dir_simu, "final_sumstats.txt")
file_stat = read.table(file_path, header = TRUE) #Lecture fichier stat résumées
#Suppression 1ère colonne, contenant uniquement chiffres pour tri en bash.
mat_stat[row_min_simu:row_max_simu,3:63] = file_stat
#Ajout de la contribution s1,0 initiale pour pouvoir calculer les delta de proportion d'admixture ultérieurement
file_path_s1.0 = str_c(dir_ne, dir_simu, str_c("simu_",nb_simu,".txt") )
file_stat_s1.0 = read.table(file_path_s1.0, header = TRUE)
mat_stat[row_min_simu:row_max_simu,64] = file_stat_s1.0$s1.0[1]
}
#print(dim(mat_stat))
}
}
#Attribution nom colonnes selon appelation stat par MetHis
colnames(mat_stat) = c("Ne", "simu", names(file_stat), "s1.0")#,"cpt")#, "delta.adm.props")
#Passage des Ne en facteur
mat_stat$Ne = factor(as.factor(mat_stat$Ne), levels = as.character(seq_ne))
#Passage simulations en entier
mat_stat$simu = as.integer(mat_stat$simu)
#Passage générations en entier
mat_stat$Generation = as.integer(mat_stat$Generation)
#Passage stats en double
for (numcol in 4:ncol(mat_stat)) {
mat_stat[,numcol] = as.double(mat_stat[,numcol])
}
# print(dim(mat_stat))
return(mat_stat)
}
Calcul du delta des moyennes de proportion d’admixture par rapport à la moyenne attendue
var.adm = function(s1, gen){
return( (s1*(1-s1))/(2**gen) )
}
delta.adm.props = function(mat){
mat$delta.mean.adm.props = NA
mat$delta.var.adm.props = NA
for (i in 1:nrow(mat)) {
ref_adm = mat$s1.0[i]
mat$delta.mean.adm.props[i] = mat$mean.adm.props[i] - ref_adm
ref_var = var.adm(s1 = ref_adm, gen = mat$Generation[i])
mat$delta.var.adm.props[i] = mat$var.adm.props[i] - ref_var
}
return(mat)
}
Dans le cas où l’on travaille sur de très nombreuses simulations correspondant à un effectif efficace, cette fonction permet de calculer les statistiques moyennées pour les valeurs de Ne étudiées.
#Passage en moyenne
data.frame.mean = function(df_stat, max_gen, seq_ne, col_mean, col_cstt){
if (length(col_mean)+length(col_cstt) != ncol(df_stat)) {
stop("all columns must be specified")
}
lrow = length(seq_ne)
df_mean = matrix(data = NA, nrow = lrow*max_gen, ncol = ncol(df_stat))
df_mean = as.data.frame(df_mean)
cpt = 0
for (ne in seq_ne) {
df_tmp = df_stat[which(df_stat$Ne == ne),]
row_min = cpt*max_gen+1
row_max = (cpt+1)*max_gen
for (cstt in col_cstt) {
# print(as.character(df_tmp[1, cstt]))
df_mean[row_min:row_max, cstt] = as.character(df_tmp[1, cstt])
}
cpt = cpt+1
for (gen in 0:(max_gen-1) ) {
df_mean[row_min+gen, 3] = gen
vec_tmp = apply(df_tmp[which(df_tmp$Generation == gen), col_mean], 2, mean)
df_mean[row_min+gen, col_mean] = vec_tmp
}
}
colnames(df_mean) = colnames(df_stat)
df_mean$Ne = as.factor(df_mean$Ne)
df_mean$Generation = as.integer(df_mean$Generation)
return(df_mean)
}
data.frame.mean.u = function(df_stat, max_gen, seq_combi, seq_u, col_mean, col_cstt){
for (u in seq_u){
if (u == seq_u[1]) {
df_tot = data.frame.mean(df_stat = df_stat[which(df_stat$U == u),],
max_gen = max_gen, seq_ne = seq_combi,
col_mean = col_mean, col_cstt = col_cstt)
}else{
df_tot = rbind(df_tot, data.frame.mean(df_stat = df_stat[which(df_stat$U == u),],
max_gen = max_gen, seq_ne = seq_combi,
col_mean, col_cstt))
}
}
df_tot$U = as.factor(as.numeric(df_tot$U))
return(df_tot)
}
data.frame.mean.bottle = function(df_stat, max_gen, seq_combi, seq_u,
seq_alpha, seq_bottle, col_mean, col_cstt){
for (alpha in seq_alpha) {
df_alpha = df_stat[which(df_stat$alpha == alpha),]
for (u in seq_u) {
df_u = df_alpha[which(df_alpha$U == u),]
for (bottle in seq_bottle) {
df_bot = df_u[which(df_u$time_botl == bottle),]
if (alpha == seq_alpha[1] && u == seq_u[1] && bottle == seq_bottle[1]) {
df_tot = data.frame.mean(df_stat = df_bot, max_gen = max_gen,
seq_ne = seq_combi, col_mean = col_mean,
col_cstt = col_cstt)
}else{
df_tot = rbind(df_tot, data.frame.mean(df_stat = df_bot,
max_gen = max_gen,
seq_ne = seq_combi,
col_mean = col_mean,
col_cstt = col_cstt))
}
}
}
}
df_tot$U = as.factor(as.numeric(df_tot$U))
df_tot$alpha = as.factor(as.numeric(df_tot$alpha))
df_tot$time_botl = as.factor(df_tot$time_botl)
df_tot$bott_u = as.factor(df_tot$bott_u)
df_tot$alpha_u = as.factor(df_tot$alpha_u)
return(df_tot)
}
Fonction d’extraction d’une sous-matrice à partir d’une matrice principale en fonction de valeurs choisies de Ne devant être présentes dans la matrice initiale.
extract_sub_mat = function(all_mat, list_ne){
all_rows = c()
for (ne in as.character(list_ne)) {
all_rows = append(all_rows, which(all_mat[,1] == ne))
}
return(all_mat[all_rows,])
}
Fonction d’affichage des graphiques en 2D, avec un background blanc et les lignes horizontales de quadrillage affichées. Cette fonction permet d’afficher les lignes reliant les points ou les lignes de régression. Il est possible de passer en paramètres la variable selon laquelle les points seront colorés, les groupes selon lesquels les points sont reliés, l’affichage ou non de la légende (l’affichage se fait par défaut), la taille des points et le type de ligne utilisé.
#Plot function
#Affichage d'une stat au cours des générations
plot_stat_gen = function(df, gen, stat, group_col = c(), color_col, titre, ligne = TRUE, legd = TRUE,
size_point = 0.1, line_t = "solid"){
p = ggplot(df, aes(x = gen, y = stat,
group = group_col, color = color_col)) + ggtitle(titre)
if (ligne) {
p = p + geom_point(size = size_point, show.legend = legd)+
geom_line(aes(linetype = line_t), show.legend = legd)
}else{
#smooth de la courbe pour obtenir régression et tempérer variabilité due au TCL
p = p + geom_point(size = size_point, show.legend = legd)+ geom_smooth(show.legend = legd)
}
p = p + theme_classic()
p = p + theme(panel.grid.major.y = element_line(size = 0.5,
linetype = 'solid',
colour = "#BABABA"),
panel.grid.minor.y = element_line(size = 0.25,
linetype = 'solid',
colour = "#CECECE"))
return(p)
}
Cette fonction d’affichage de graphique reprend la précédente, en affichant uniquement les courbes de régression sans les points.
#Plot function
#Affichage d'une stat au cours des générations
plot_without_point = function(df, gen, stat, color_col, titre, legd = TRUE){
p = ggplot(df, aes(x = gen, y = stat, color = color_col)) + ggtitle(titre)
#smooth de la courbe pour obtenir régression et tempérer variabilité due au TCL
p = p + geom_smooth(show.legend = legd)
p = p + theme_classic()
p = p + theme(panel.grid.major.y = element_line(size = 0.5,
linetype = 'solid',
colour = "#BABABA"),
panel.grid.minor.y = element_line(size = 0.25,
linetype = 'solid',
colour = "#CECECE"))
#print(p)
return(p)
}
Cette fonction permet de rajouter à un plot existant : - des lignes verticales (pour les goulots d’étranglement par exemple) - un titre aux légendes de couleurs et de types de ligne - de modifier les types de lignes utilisés
improve_plot_bottle = function(p, vec_abline, lab_col, lab_line, vec_linetype){
for (rg_abl in 1:length(vec_abline)) {
p = p + geom_vline(xintercept = vec_abline[rg_abl], size = 0.4, color = "orange", linetype = "solid")
}
p = p + labs(color = lab_col, linetype = lab_line) + scale_linetype_manual(values=vec_linetype)
return(p)
}
Fonction pour une représentation graphique à l’aide d’un gif, pour suivre l’évolution de la distribution des proportions d’admixture dans la population.
hist.gif.props.adm = function(mat, select_stat, select_seq, title, title_leg,
seq_gen = c(seq(1,10,1), seq(11,101,10)),
vec_col = c("#A5260A", "#F6DC12", "#87E990", "#FFE4C4",
"#E67E30", "#79F8F8", "#A10684", "#8B6C42")){
png(file="example%03d.png", width=600, heigh=400)
for (line in seq_gen) {
cpt_col = 1
for (var in select_seq) {
rg_select = which(colnames(mat) == select_stat)
mat_tmp = mat[which(mat[,rg_select] == var),]
vec_y = c(mat_tmp$perc0.adm.props[line],
mat_tmp$perc10.adm.props[line],
mat_tmp$perc20.adm.props[line],
mat_tmp$perc30.adm.props[line],
mat_tmp$perc40.adm.props[line],
mat_tmp$perc50.adm.props[line],
mat_tmp$perc60.adm.props[line],
mat_tmp$perc70.adm.props[line],
mat_tmp$perc80.adm.props[line],
mat_tmp$perc90.adm.props[line],
mat_tmp$perc100.adm.props[line])
if (var == select_seq[1]) {
hist(vec_y, breaks = vec_y, main = line-1, col = vec_col[cpt_col], xlim = c(-0.2,1.2))
cpt_col = cpt_col+1
}else{
hist(vec_y, breaks = vec_y, add = TRUE, col = vec_col[cpt_col])
cpt_col = cpt_col+1
}
}
legend("topleft", title=title_leg,legend = select_seq,
fill = vec_col[1:cpt_col], cex=0.8)
}
dev.off()
system(str_c("convert -delay 120 *.png ",title,".gif"))
file.remove(list.files(pattern=".png"))
}
print(difftime(Sys.time(), time.ini))
## Time difference of 0.121845 secs
Calcul des matrices de stats pour de populations constantes de tailles initiales différentes.
seq_1024_16384 = 2**seq(10,14,1)
seq_50_1000 = seq(50,1000,1)
seq_tot = sort(c(seq_50_1000, seq_1024_16384))
# mat tot cstt
mat_tot = data.frame.stat(seq_tot, 101, 1)
# Ne de 50 à 200 par pas de 10 puis de 200 à 1000 par pas de 100 (mat1 cstt)
mat_50_1000 = extract_sub_mat(mat_tot, c(seq(50,200,10), seq(200,1000,100)))
# Ne de 1024 à 16384 avec un pas de x2 (mat2 cstt)
mat_1024_16384 = extract_sub_mat(mat_tot, 2**seq(6,14,1))
# Ne de 64 à 16384 avec un pas de x2 + Ne de 100 et de 1000 (mat3 cstt)
mat_64_100_16384 = extract_sub_mat(mat_tot, c(2**seq(6,14,1), 100, 1000) )
Graphiques obtenus : affichage de l’évolution d’une statistique donnée en fonction des générations colorée selon la Ne initiale.
# Affichage de l'hétérozygotie pour la matrice totale (mat tot cstt)
p = plot_stat_gen(mat_tot, mat_tot$Generation,
mat_tot$mean.het.adm, mat_tot$Ne, mat_tot$Ne,ligne = T,
"Moyenne He en fonction des générations\n selon différentes Ne initiales",
legd = FALSE)
p = p+geom_hline(yintercept = mat_50_1000$mean.het.s1[1], color = "red")
p = p+geom_hline(yintercept = mat_50_1000$mean.het.s2[1], color = "black")
print(p)
# Affichage de l'hétérozygotie pour mat1 cstt
p = plot_stat_gen(mat_50_1000, mat_50_1000$Generation,
mat_50_1000$mean.het.adm, mat_50_1000$Ne, mat_50_1000$Ne,ligne = T,
"Moyenne He en fonction des générations\nselon différentes Ne initiales",
legd = TRUE)
p = p+geom_hline(yintercept = mat_50_1000$mean.het.s1[1], color = "red")
p = p+geom_hline(yintercept = mat_50_1000$mean.het.s2[1], color = "black")
print(p)
plot_stat_gen(mat_tot, mat_tot$Generation,
mat_tot$Fst.s1.adm, mat_tot$Ne,mat_tot$Ne,
"Fst s1 adm en fonction des générations\n selon différentes Ne initiales",
legd = FALSE)
plot_stat_gen(mat_50_1000, mat_50_1000$Generation,
mat_50_1000$Fst.s1.adm, mat_50_1000$Ne,mat_50_1000$Ne,
"Fst s1 adm en fonction des générations\n selon différentes Ne initiales")
plot_stat_gen(mat_1024_16384, mat_1024_16384$Generation,
mat_1024_16384$Fst.s1.adm, mat_1024_16384$Ne,mat_1024_16384$Ne,ligne = T,
"Fst s1 adm en fonction des générations\n selon différentes Ne initiales")
Fonction d’écriture de plot dans des fichiers, en choisissant la statistique enregistrée.
write_cstt_plot = function(name_stat, mat, name_file,min_y, max_y, name_x = "Generation"){
num_col_y = which(colnames(mat) == name_stat)
num_col_x = which(colnames(mat) == name_x)
p = plot_stat_gen(mat, mat[,num_col_x],
mat[,num_col_y], mat$Ne, mat$Ne,ligne = T,
str_c(name_stat, " en fonction des générations\n selon différentes Ne initiales"))
p = p+ylim(min_y, max_y)+ labs(color = "Ne") + guides(linetype = FALSE)
p = p+xlab(name_x)+ylab(name_stat)
png(filename = str_c("../../../Images/", name_file, ".png"), width = 1000, height = 800)
print(p)
dev.off()
ggsave(filename = str_c("../../../Images/", name_file, ".eps"), plot = p)
}
Ajout delta mean et var des proportions d’admixture pour les effectifs efficaces constants.
mat_64_100_16384 = delta.adm.props(mat_64_100_16384)
Écriture de différents graphiques pour différentes statistiques choisies. La matrice choisie est la mat3 cstt pour des raisons de lisibilité. Cette matrice permet également d’observer l’évolution des statistiques choisies pour des Ne = 100 et Ne = 1000, valeurs qui seront reprises dans d’autres scénarios ultérieurement.^
write_cstt_plot("mean.het.adm", mat_64_100_16384, "moyenne_heterozygotie/moy_He_ne_cst", 0.01, 0.062)
## Saving 7 x 5 in image
write_cstt_plot("Fst.s1.adm", mat_64_100_16384, "Fsts1adm/fst_s1_adm_ne_cst", -0.01, 0.5)
## Saving 7 x 5 in image
write_cstt_plot("Fst.s2.adm", mat_64_100_16384, "Fsts2adm/fst_s2_adm_ne_cst", -0.01, 0.5)
## Saving 7 x 5 in image
write_cstt_plot("mean.F.adm", mat_64_100_16384, "Fadm/f_adm_ne_cst", -0.03, 0.03)
## Warning: Removed 6 rows containing missing values (geom_point).
## Warning: Removed 1 row(s) containing missing values (geom_path).
## Saving 7 x 5 in image
## Warning: Removed 6 rows containing missing values (geom_point).
## Warning: Removed 1 row(s) containing missing values (geom_path).
write_cstt_plot("delta.mean.adm.props", mat_64_100_16384, "delta_mean_adm/delta_mean_ne_cstt", -0.25, 0.25)
## Saving 7 x 5 in image
write_cstt_plot("delta.var.adm.props", mat_64_100_16384, "delta_var_adm/delta_var_ne_cstt", -0.2, 0.05)
## Saving 7 x 5 in image
setwd("../../../Images/")
hist.gif.props.adm(mat = mat_64_100_16384, select_stat = "Ne", select_seq = c(64,100,256,2048),
title = "adm.props/hist_cstt", title_leg = "Ne")
## [1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [16] TRUE TRUE TRUE TRUE TRUE
print(difftime(Sys.time(), time.ini))
## Time difference of 2.236184 mins
Réutilisation de la fonction de lecture de stats résumées pour des populations croissantes. Celles-ci sont enregistrées dans des répertoires de la forme “Ne100-XXX/Ne100-1000/simu_1” Fonction de lecture des stats résumées. La fonction est prévue pour lire les fichiers de stats résumées dans des répertoires de la forme “NeXXX/simu_XXX” (avec XXX remplacés par les valeurs de Ne et de simulations). Elle renvoie les résultats sous forme de data frame en ajoutant aux statistiques résumées les colonnes correspondant à l’effectif efficace et la simulation.
data.frame.ne.gen = function(seq_ne){
for (ne in seq_ne) {
if (dir.exists(str_c("Ne", ne))) {
dir_ne = str_c("Ne", ne, "/simu_1/")
setwd(dir_ne)
file_tmp = read.table("simu_1.par", header = TRUE)
if (ne == seq_ne[1]) {
mat_ne_inc = as.data.frame(file_tmp[,2])
}else{
mat_ne_inc = cbind(mat_ne_inc, file_tmp[,2])
}
setwd("../../")
}
}
return(mat_ne_inc)
}
data.frame.increase = function(seq_ne_ini, seq_combi, max_simu = 1){
mat_inc = c()
for (ne_ini in seq_ne_ini) {
if (dir.exists(str_c("Ne", ne_ini, "-XXX/"))) {
dir_ne = str_c("Ne", ne_ini, "-XXX/")
setwd(dir_ne)
motif_detect = str_c("^",ne_ini,"-")
vec_ne_tmp = seq_combi[which(str_detect(seq_combi, motif_detect))]
mat_tmp = data.frame.stat(seq_ne = vec_ne_tmp, max_gen = 101, max_simu = max_simu)
mat_ne = data.frame.ne.gen(seq_ne = vec_ne_tmp)
# print(mat_ne)
vec_tmp = as.integer(unlist(str_split(mat_tmp$Ne, "-")))
mat_tmp$ne_gen = unlist(mat_ne, use.names = F)
# print(length(vec_tmp))
mat_tmp$n0 = vec_tmp[seq(1,length(vec_tmp),2)]
mat_tmp$nf = vec_tmp[seq(2,length(vec_tmp),2)]
if (ne_ini == seq_ne_ini[1]) {
mat_inc = mat_tmp
}else{
mat_inc = rbind(mat_inc, mat_tmp)
}
setwd("../")
#print(dim(mat_pop_inc))
}
}
return(na.omit(mat_inc))
}
# mat_pop_nu = data.frame(mat_pop_nu,
# ne_gen = unlist(mat_ne_inc , use.names = F),
# n0 = rep(NA, nrow(mat_pop_nu)),
# nf = rep(NA, nrow(mat_pop_nu)) )
#
# for (i in 1:nrow(mat_pop_nu)) {
# vec_tmp = as.integer(unlist(str_split(mat_pop_nu$Ne, "-")))
# mat_pop_nu$n0[i] = vec_tmp[1]
# mat_pop_nu$nf[i] = vec_tmp[2]
# }
#
# mat_ne_inc = data.frame()
# cpt = 1
Lecture des data frame pour les combinaison de Ne entre ne0 (ne_ini) et nef (ne_fin). On extrait ensuite les stats pour lesquelles on a un nef de 10000.
setwd("../new_methis_pop_increase_50000_snp/")
seq_ne_ini = seq(50,100,2)
seq_ne_fin = c(seq(100, 10000, 200), 10000)
seq_combi = as.data.frame(outer(seq_ne_ini, seq_ne_fin, FUN="paste", sep="-"))
seq_combi = as.data.frame.vector(unlist(seq_combi))[[1]]
mat_pop_inc = data.frame.increase(seq_ne_ini, seq_combi)
mat_end_10000 = extract_sub_mat(mat_pop_inc,
as.data.frame(outer(seq(50,100,4),
c("10000"),
FUN="paste",
sep="-"))[,1])
plot_stat_gen(df = mat_pop_inc, gen = mat_pop_inc$Generation,
stat = mat_pop_inc$f3, group_col = mat_pop_inc$Ne,color_col = mat_pop_inc$Ne,
titre = "Stat en fonction des générations\n selon différentes Ne initiales",
legd = FALSE)
plot_stat_gen(mat_pop_inc, mat_pop_inc$Generation,
mat_pop_inc$mean.het.adm, mat_pop_inc$Ne,mat_pop_inc$Ne,
"Stat en fonction des générations\n selon différentes Ne initiales",
legd = FALSE)
plot_stat_gen(mat_end_10000, mat_end_10000$Generation,
mat_end_10000$mean.het.adm, mat_end_10000$Ne,mat_end_10000$Ne,
"Moyenne hétérozygotie en fonction des générations\n selon différentes Ne initiales",
legd = TRUE)
plot_stat_gen(mat_end_10000, mat_end_10000$Generation,
mat_end_10000$f3, mat_end_10000$Ne,mat_end_10000$Ne,
"F3 en fonction des générations\n selon différentes Ne initiales",
legd = TRUE)
print(difftime(Sys.time(), time.ini))
## Time difference of 2.886243 mins
Le U déterminant l’inflexion de la courbe de croissance est ici connu et fixé avant la simulation.
data.frame.increase.nu = function(seq_ne_ini, seq_combi, seq_nu, max_simu = 1){
for (nu in seq_nu) {
if(dir.exists(str_c("Nu", nu ))){
dir_nu = str_c("Nu", nu, "/")
setwd(dir_nu)
mat_tmp = data.frame.increase(seq_ne_ini, seq_combi, max_simu)
if (nu == seq_nu[1]) {
mat_pop = data.frame(mat_tmp, U = nu)
}else{
mat_tmp = data.frame(mat_tmp, U = nu)
mat_pop = rbind(mat_pop, mat_tmp)
}
setwd("../")
}
}
mat_pop$U = as.factor(mat_pop$U)
return(mat_pop)
}
setwd("../new_methis_pop_increase_50000_snp_Nu_var/")
#Détermination des paramètres initiaux
seq_nu = sprintf("%.2f", seq(0.05, 0.45, 0.05))
seq_ne_ini = seq(50,100,10)
seq_ne_fin = seq(100, 10000, 200)
seq_combi = as.data.frame(outer(seq_ne_ini, seq_ne_fin, FUN="paste", sep="-"))
seq_combi = as.data.frame.vector(unlist(seq_combi))[[1]]
mat_pop_nu = data.frame.increase.nu(seq_ne_ini = seq_ne_ini, seq_combi = seq_combi, seq_nu = seq_nu)
plot_without_point(mat_pop_nu, mat_pop_nu$Generation,
mat_pop_nu$mean.het.adm, mat_pop_nu$Ne,
"Moyenne hétérozygotie en fonction des générations\n selon différentes Ne initiales\n et couleurs selon U",
legd = FALSE)
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
plot_without_point(mat_pop_nu, mat_pop_nu$Generation,
mat_pop_nu$mean.het.adm, mat_pop_nu$U,
"Moyenne hétérozygotie en fonction des générations\n selon différentes Ne initiales\n et couleurs selon nu",
legd = TRUE)
## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
# plot_stat_gen(mat_pop_nu, mat_pop_nu$Generation,
# mat_pop_nu$mean.het.adm, mat_pop_nu$Ne,
# "Stat en fonction des générations\n selon différentes Ne initiales",
# TRUE, legd = FALSE)
# plot_stat_gen(mat_pop_nu, mat_pop_nu$Generation,
# mat_pop_nu$mean.het.adm, mat_pop_nu$Nu,
# "Stat en fonction des générations\n selon différentes Ne initiales",
# TRUE, legd = FALSE)
Affichage de l’hétérozygotie pour les différentes valeurs de Nu avec N0 entre 50 et 100 (pas de 10) et Nf entre 100 et 10 000 (pas de 200).
for (nu in seq_nu) {
mat_pop_nu_tmp = mat_pop_nu[which(mat_pop_nu$U == nu),]
p = plot_stat_gen(mat_pop_nu_tmp, mat_pop_nu_tmp$Generation,
mat_pop_nu_tmp$mean.het.adm, mat_pop_nu_tmp$Ne,
mat_pop_nu_tmp$Ne,
str_c("U = ",nu), TRUE, legd = FALSE)
p = p+geom_hline(yintercept = mat_pop_nu_tmp$mean.het.s1[1],
color = "red")
p = p+geom_hline(yintercept = mat_pop_nu_tmp$mean.het.s2[1],
color = "black")
print(p)
}
remove(mat_pop_nu_tmp)
Affichage en surface de l’hétérozygotie à la génération 100 en fonction de Ne0 et Nef selon différents U.
for (nu in seq_nu) {
mat_pop_nu_tmp = mat_pop_nu[which(mat_pop_nu$U == nu),]
list_ne = mat_pop_nu_tmp$Ne[which(mat_pop_nu_tmp$Generation == 100)]
list_ne = as.integer(unlist(str_split(list_ne, "-")))
list_ne0 = list_ne[seq(1, length(list_ne), 2)]
list_nef = list_ne[seq(2, length(list_ne), 2)]
mat_pop_nu_tmp = mat_pop_nu_tmp[which(mat_pop_nu_tmp$Generation == 100),]
print(wireframe(mat_pop_nu_tmp$mean.het.adm ~ list_ne0 * list_nef, data = mat_pop_nu_tmp,
drape = TRUE, zlab = "He", xlab = "Ne0", ylab = "Nef",
screen = list(z = 40, x = -60),
scales = list(z.ticks=5, arrows=F, col="black", font=1, tck=0.8),
main = str_c("He moyenne à la génération 100 en fonction de Ne0 et Nef\n U = ", nu) ))
}
remove(mat_pop_nu_tmp)
for (rg_nu in seq_nu) {
#rg_nu
df_tmp = mat_pop_nu[which(mat_pop_nu$U == rg_nu),]
scatter3D(x = df_tmp$Generation, y = df_tmp$ne_gen, z = df_tmp$mean.het.adm,
bty = "b2", pch = 16, ticktype = "detailed", main = rg_nu,
theta = 210, phi = 20, type = "s",cex = 0.2, cex.lab = 0.6, cex.axis = 0.5,
xlab = "Gen", ylab = "Ne", zlab = "Stat")
# print(wireframe(mean.het.adm ~ Generation * ne_gen,
# data = df_tmp,
# drape = TRUE, zlab = "Htz", xlab = "generation", ylab = "Ne gen",
# screen = list(z = 40, x = -60),
# scales = list(z.ticks=5, arrows=F, col="black", font=1, tck=0.8),
# main = rg_nu))
}
rm(df_tmp)
scatter3D(x = mat_pop_nu$Generation, y = mat_pop_nu$ne_gen,
z = as.numeric(as.character(mat_pop_nu$U)),
bty = "b2", pch = 16, ticktype = "detailed",
theta = 210, phi = 50, type = "s", cex.lab = 0.6, cex.axis = 0.5,
xlab = "Gen", ylab = "Ne", zlab = "Nu")
setwd("../new_methis_pop_increase_50000_snp_Nu_var/")
seq_nu_bis = sprintf("%.2f", c(0.01, 0.25, 0.49))
seq_ne_ini = c(100)
seq_ne_fin = c(200, 400, seq(1000, 10000, 1000))
seq_combi = as.data.frame(outer(seq_ne_ini, seq_ne_fin, FUN="paste", sep="-"))
seq_combi = as.data.frame.vector(unlist(seq_combi))[[1]]
mat_pop_nu_reduce = data.frame.increase.nu(seq_ne_ini, seq_combi, seq_nu_bis)
mat_pop_nu_reduce_gen_100 = mat_pop_nu_reduce[which(mat_pop_nu_reduce$Generation == 100),]
mat_pop_nu_reduce$log_mean = log(mat_pop_nu_reduce$mean.het.adm)
scatter3D(x = mat_pop_nu_reduce$Generation, y = mat_pop_nu_reduce$nf,
z = mat_pop_nu_reduce$mean.het.adm, colvar = as.numeric(mat_pop_nu_reduce$U),
col = c("#1B9E77", "#D95F02", "#7570B3"),
colkey = list(at = c(1.3, 2, 2.7),
addlines = TRUE, length = 1, width = 0.5,
labels = levels(as.factor(mat_pop_nu_reduce$U))),
ticktype = "detailed",bty = "g",pch = 16,
theta = 115, phi = 0, type = "s", cex = 0.5, cex.lab = 0.6, cex.axis = 0.5,
xlab = "generation", ylab = "Nf", zlab = "stat", clab = "U")
# scatter3D(x = df_tmp$Generation, y = vec_tot, z = df_tmp$mean.het.adm,
# bty = "b2", pch = 16, ticktype = "detailed", main = seq_nu[rg_nu],
# theta = 210, phi = 20, type = "s", cex.lab = 0.6, cex.axis = 0.5,
# xlab = "Gen", ylab = "Ne", zlab = "Stat")
setwd("../new_methis_pop_increase_50000_snp_Nu_var/")
seq_ne_ini3 = seq(100,100,0)
seq_multi3 = c(2,4,10,100)
seq_ne_fin3 = seq_multi3*seq_ne_ini3
seq_combi3 = as.data.frame(outer(seq_ne_ini3, seq_ne_fin3, FUN="paste", sep="-"))
seq_combi3 = as.data.frame.vector(unlist(seq_combi3))[[1]]
mat_pop_nu_reduce_3 = data.frame.increase.nu(seq_ne_ini3, seq_combi3, seq_nu_bis)
write_increase_plot = function(name_stat, mat, name_file, min_y, max_y, name_x = "Generation"){
num_col_y = which(colnames(mat) == name_stat)
num_col_x = which(colnames(mat) == name_x)
p = plot_stat_gen(df = mat, gen = mat[,num_col_x],
stat = mat[,num_col_y], group_col = c(),
color_col = mat$Ne, line_t = mat$U,
titre = str_c(name_stat," en fonction des générations\n selon différentes croissances de populations"),
legd = TRUE)
p = improve_plot_bottle(p =p, vec_abline = c(), lab_col = "ne0-nef",
lab_line = "U", vec_linetype = c("dotted", "longdash", "solid"))
p = p + ylim(min_y, max_y)
p = p+xlab(name_x)+ylab(name_stat)
png(filename = str_c("../../../Images/", name_file,".png"), width = 1000, height = 800)
print(p)
dev.off()
ggsave(filename = str_c("../../../Images/", name_file, ".eps"), plot = p)
}
Ajout delta mean et var des proportions d’admixture pour les populations croissantes.
mat_pop_nu_reduce_3 = delta.adm.props(mat_pop_nu_reduce_3)
write_increase_plot("mean.het.adm", mat_pop_nu_reduce_3, "moyenne_heterozygotie/moy_He_pop_inc_nu_var",0.01,0.062)
## Saving 7 x 5 in image
write_increase_plot("Fst.s1.adm", mat_pop_nu_reduce_3, "Fsts1adm/fst_s1_adm_pop_inc_nu_var", 0, 0.5)
## Warning: Removed 2 rows containing missing values (geom_point).
## Warning: Removed 2 row(s) containing missing values (geom_path).
## Saving 7 x 5 in image
## Warning: Removed 2 rows containing missing values (geom_point).
## Warning: Removed 2 row(s) containing missing values (geom_path).
write_increase_plot("Fst.s2.adm", mat_pop_nu_reduce_3, "Fsts2adm/fst_s2_adm_pop_inc_nu_var", 0, 0.5)
## Saving 7 x 5 in image
write_increase_plot("mean.F.adm", mat_pop_nu_reduce_3, "Fadm/f_adm_pop_inc_nu_var", -0.03, 0.03)
## Warning: Removed 1 rows containing missing values (geom_point).
## Saving 7 x 5 in image
## Warning: Removed 1 rows containing missing values (geom_point).
write_increase_plot("delta.mean.adm.props", mat_pop_nu_reduce_3, "delta_mean_adm/delta_mean_pop_inc_nu_va", -0.25, 0.25)
## Saving 7 x 5 in image
write_increase_plot("delta.var.adm.props", mat_pop_nu_reduce_3, "delta_var_adm/delta_var_pop_inc_nu_va", -0.2, 0.05)
## Saving 7 x 5 in image
setwd("../../../Images/")
mat_pop_nu_reduce_100_10000 = mat_pop_nu_reduce_3[which(mat_pop_nu_reduce_3$Ne == "100-10000"),]
hist.gif.props.adm(mat = mat_pop_nu_reduce_100_10000, select_stat = "U", select_seq = seq_nu_bis,
title = "adm.props/hist_inc", title_leg = "U")
## logical(0)
print(difftime(Sys.time(), time.ini))
## Time difference of 4.53561 mins
data.frame.bottleneck = function(lst_ne0, lst_nef, lst_alpha, lst_nu, lst_bottle, max_simu = 1){
lst_combi = as.data.frame(outer(lst_ne0, lst_nef, FUN="paste", sep="-"))
lst_combi = as.data.frame.vector(unlist(lst_combi))[[1]]
cpt = 0
for (alpha in lst_alpha) {
for (nu in lst_nu) {
for (bottle in lst_bottle) {
change_dir = str_c("alpha", alpha, "/Nu", nu, "/bottle",bottle,"/")
# print(change_dir)
if(dir.exists(change_dir)){
# print(getwd())
setwd(change_dir)
mat_tmp = data.frame.increase(lst_ne0, lst_combi, max_simu)
if (!is.null(nrow(mat_tmp))) {
col_alpha = rep(alpha, nrow(mat_tmp))
col_nu = rep(nu, nrow(mat_tmp))
col_bottle = rep(bottle, nrow(mat_tmp))
col_bott_u = rep(str_c(bottle, "/", nu), nrow(mat_tmp))
# col_bott_u = as.numeric(rep(bottle, nrow(mat_tmp)))/as.numeric(rep(nu, nrow(mat_tmp)))
# col_bott_u = ceiling(col_bott_u/1)*1
col_alpha_u = rep(str_c(nu, "/", alpha), nrow(mat_tmp))
# col_alpha_u = as.numeric(rep(nu, nrow(mat_tmp)))/as.numeric(rep(alpha, nrow(mat_tmp)))
# col_alpha_u = ceiling(col_alpha_u/0.001)*0.001
mat_tmp = data.frame(mat_tmp, alpha = col_alpha,
U = col_nu, time_botl = col_bottle,
bott_u = col_bott_u, alpha_u = col_alpha_u)
cpt = cpt+1
}
# print(dim(mat_tmp))
if (alpha==lst_alpha[1] & nu==lst_nu[1] & bottle==lst_bottle[1]) {
mat_bot = mat_tmp
}else{
mat_bot = rbind(mat_bot, mat_tmp)
}
setwd("../../../")
}
}
}
}
if(!is.null(mat_bot)){
mat_bot$alpha = as.factor(mat_bot$alpha)
mat_bot$U = as.factor(mat_bot$U)
mat_bot$time_botl = as.factor(mat_bot$time_botl)
mat_bot$bott_u = as.factor(mat_bot$bott_u)
mat_bot$alpha_u = as.factor(mat_bot$alpha_u)
}
return(mat_bot)
}
setwd("../new_methis_bottleneck_50000_snp/")
mat_bottle = data.frame.bottleneck(lst_ne0 = c(1000), lst_nef = c(1000),
lst_alpha = c(0.1, 0.5, 0.9),
lst_nu = c(0.01, 0.25, 0.49),
lst_bottle = seq(10, 90, 10))
mat_bottle_time_50 = mat_bottle[which(mat_bottle$time_botl == 50),]
mat_bottle_time_20 = mat_bottle[which(mat_bottle$time_botl == 20),]
mat_bottle_time_80 = mat_bottle[which(mat_bottle$time_botl == 80),]
boucle_plot_bottle = function(mat, vec_stat, name_stat, vec_bott,size_pop = 1000){
for (bott in vec_bott) {
new_mat = mat[which(mat$time_botl == bott),]
new_vec = vec_stat[which(mat$time_botl == bott)]
p = plot_stat_gen(new_mat, new_mat$Generation,
new_vec, c(),
new_mat$alpha, size_point = 0,
str_c(name_stat, " en fonction des générations\nbottleneck ",
bott, " N0 ",size_pop," Nf ",size_pop),
TRUE,legd = TRUE, line_t = new_mat$U)
p = improve_plot_bottle(p = p, vec_abline = c(1,bott+1), lab_col = "alpha",
lab_line = "U", vec_linetype = c("solid", "longdash", "dotted"))
print(p)
}
}
plot_without_point(mat_bottle, mat_bottle$Generation,
mat_bottle$mean.het.adm, mat_bottle$alpha,
"Stat en fonction des générations\n selon différents bottleneck",
legd = TRUE)
## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
Moyenne de l’hétérozygotie
boucle_plot_bottle(mat_bottle, mat_bottle$mean.het.adm, "Heterozygotie", c(20,50,80))
Fst
boucle_plot_bottle(mat_bottle, mat_bottle$Fst.s1.adm, "Fst s1 adm", c(20,50,80))
boucle_plot_bottle(mat_bottle, mat_bottle$Fst.s2.adm, "Fst s2 adm", c(20,50,80))
boucle_plot_bottle(mat_bottle, mat_bottle$mean.F.adm, "F mean", c(20,50,80))
mini_mat_bottle = rbind(mat_bottle_time_20, mat_bottle_time_50, mat_bottle_time_80)
p = plot_stat_gen(mat_bottle, mat_bottle$Generation,
mat_bottle$mean.het.adm, c(),
mat_bottle$bott_u, line_t = mat_bottle$alpha,
"Stat en fonction des générations\n selon différents bottleneck",TRUE,
legd = TRUE)
p = p+ theme(legend.position="bottom", legend.box = "horizontal")
print(p)
p = plot_stat_gen(mini_mat_bottle, mini_mat_bottle$Generation,
mini_mat_bottle$mean.het.adm, c(),
mini_mat_bottle$bott_u, line_t = mini_mat_bottle$alpha,
"Stat en fonction des générations\n selon différents bottleneck",TRUE,
legd = TRUE)
improve_plot_bottle(p =p, vec_abline = c(21, 51, 81), lab_col = "tpsbott/U",
lab_line = "alpha", vec_linetype = c("solid", "longdash", "dotted"))
write_bottle_plot = function(name_stat, mat, name_file, min_y, max_y, name_x = "Generation"){
num_col_y = which(colnames(mat) == name_stat)
num_col_x = which(colnames(mat) == name_x)
p = plot_stat_gen(mat, mat[,num_col_x],
mat[,num_col_y], c(),
mat$alpha_u, line_t = mat$time_botl,
str_c(name_stat, " en fonction des générations\n selon différents bottleneck\n (Ne0 = Nef = 1000)"),TRUE,
legd = TRUE)
p = improve_plot_bottle(p =p, vec_abline = c(21, 51, 81), lab_col = "U/alpha",
lab_line = "tpsbott", vec_linetype = c("dotted", "longdash", "solid"))
p = p+ylim(min_y, max_y)
p = p+xlab(name_x)+ylab(name_stat)
png(filename = str_c("../../../Images/",name_file,".png"), width = 1000, height = 800)
print(p)
dev.off()
ggsave(filename = str_c("../../../Images/", name_file, ".eps"), plot = p)
}
p = plot_stat_gen(mini_mat_bottle, mini_mat_bottle$Generation,
mini_mat_bottle$mean.het.adm, c(),
mini_mat_bottle$alpha_u, line_t = mini_mat_bottle$time_botl,
"Moyenne He en fonction des générations\n selon différents bottleneck\n (Ne0 = Nef = 1000)",
TRUE, legd = TRUE)
p = improve_plot_bottle(p =p, vec_abline = c(21, 51, 81), lab_col = "U/alpha",
lab_line = "tpsbott", vec_linetype = c("dotted", "longdash", "solid"))
print(p)
p = plot_stat_gen(mini_mat_bottle, mini_mat_bottle$Generation,
mini_mat_bottle$Fst.s1.adm, c(),
mini_mat_bottle$alpha_u, line_t = mini_mat_bottle$time_botl,
"Fst s1 adm en fonction des générations\n selon différents bottleneck\n (Ne0 = Nef = 1000)",
TRUE, legd = TRUE)
p = improve_plot_bottle(p =p, vec_abline = c(21, 51, 81), lab_col = "U/alpha",
lab_line = "tpsbott", vec_linetype = c("dotted", "longdash", "solid"))
print(p)
p = plot_stat_gen(mini_mat_bottle, mini_mat_bottle$Generation,
mini_mat_bottle$mean.F.adm, c(),
mini_mat_bottle$alpha_u, line_t = mini_mat_bottle$time_botl,
"Moyenne F en fonction des générations\n selon différents bottleneck\n (Ne0 = Nef = 1000)",
TRUE, legd = TRUE)
p = improve_plot_bottle(p =p, vec_abline = c(21, 51, 81), lab_col = "U/alpha",
lab_line = "tpsbott", vec_linetype = c("dotted", "longdash", "solid"))
print(p)
# p = plot_stat_gen(mini_mat_bottle, mini_mat_bottle$Generation,
# mini_mat_bottle$delta.adm.props, c(),
# mini_mat_bottle$alpha_u, line_t = mini_mat_bottle$time_botl,
# "delta props en fonction des générations\n selon différents bottleneck\n (Ne0 = Nef = 1000)",TRUE,
# legd = TRUE)
# p = improve_plot_bottle(p =p, vec_abline = c(21, 51, 81), lab_col = "U/alpha",
# lab_line = "tpsbott", vec_linetype = c("dotted", "longdash", "solid"))
# print(p)
Ajout delta mean et var des proportions d’admixture pour les effectifs efficaces constants.
mini_mat_bottle = delta.adm.props(mini_mat_bottle)
write_bottle_plot("mean.het.adm", mini_mat_bottle, "moyenne_heterozygotie/moy_He_bottleneck",0.01,0.062)
## Saving 7 x 5 in image
write_bottle_plot("Fst.s1.adm", mini_mat_bottle, "Fsts1adm/fst_s1_adm_bottleneck", 0, 0.5)
## Warning: Removed 1 rows containing missing values (geom_point).
## Warning: Removed 1 row(s) containing missing values (geom_path).
## Saving 7 x 5 in image
## Warning: Removed 1 rows containing missing values (geom_point).
## Warning: Removed 1 row(s) containing missing values (geom_path).
write_bottle_plot("Fst.s2.adm", mini_mat_bottle, "Fsts2adm/fst_s2_adm_bottleneck", 0, 0.5)
## Saving 7 x 5 in image
write_bottle_plot("mean.F.adm", mini_mat_bottle, "Fadm/f_adm_bottleneck", -0.03, 0.03)
## Saving 7 x 5 in image
write_bottle_plot("delta.mean.adm.props", mini_mat_bottle, "delta_mean_adm/delta_mean_bottleneck", -0.25, 0.25)
## Saving 7 x 5 in image
write_bottle_plot("delta.var.adm.props", mini_mat_bottle, "delta_var_adm/delta_var_bottleneck", -0.2, 0.05)
## Saving 7 x 5 in image
mini_mat_bottle$log.delta.var = log2(abs(mini_mat_bottle$delta.var.adm.props))
write_bottle_plot("log.delta.var", mini_mat_bottle, "delta_var_adm/log_delta_var_bottleneck", -15, 0)
## Warning: Removed 8 rows containing missing values (geom_point).
## Saving 7 x 5 in image
## Warning: Removed 8 rows containing missing values (geom_point).
vec_tmp = c()
vec_tmp[which(mat_bottle$alpha == 0.1)] = 2
vec_tmp[which(mat_bottle$alpha == 0.5)] = 3
vec_tmp[which(mat_bottle$alpha == 0.9)] = 4
scatter3D(x = mat_bottle$Generation, y = as.numeric(as.character(mat_bottle$time_botl)),
z = mat_bottle$mean.het.adm, colvar = as.numeric(vec_tmp),
ticktype = "detailed",bty = "b2",pch = 16,
theta = 70, phi = 0, type = "s",cex = 0.3, cex.lab = 0.6, cex.axis = 0.5,
xlab = "generation", ylab = "bottleneck", zlab = "He moyen",
col = c("#1B9E77", "#D95F02", "#7570B3"),
colkey = list(at = c(2.3, 3, 3.7),
addlines = TRUE, length = 1, width = 0.5,
labels = c("0.1", "0.5", "0.9")),
clab = "alpha",
main = "He moyenne en fonction des générations et\n du temps de bottleneck")
remove(vec_tmp)
setwd("../../../Images/")
seq_alpha = c(0.1, 0.5, 0.9)
mat_bottle_time_80_u_0.01 = mat_bottle_time_80[which(mat_bottle_time_80$U == 0.01),]
hist.gif.props.adm(mat = mat_bottle_time_80_u_0.01, select_stat = "alpha", select_seq = seq_alpha,
title = "adm.props/hist_bottle", title_leg = "alpha")
## [1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [16] TRUE TRUE TRUE TRUE TRUE
mat_bottle_time_80_u_0.01 = delta.adm.props(mat_bottle_time_80_u_0.01)
# plot(mat_bottle_time_80_u_0.01$delta.mean.adm.props~mat_bottle_time_80_u_0.01$Generation, type = "l")
# plot(mat_bottle_time_80_u_0.01$delta.var.adm.props~mat_bottle_time_80_u_0.01$Generation, type = "l")
plot_stat_gen(mat_bottle_time_80_u_0.01, mat_bottle_time_80_u_0.01$Generation,
mat_bottle_time_80_u_0.01$delta.mean.adm.props,
mat_bottle_time_80_u_0.01$alpha, mat_bottle_time_80_u_0.01$alpha,ligne = T,
"Delta mean adm props en fonction des générations\n bottleneck = 80, U = 0.01",
legd = TRUE)
plot_stat_gen(mat_bottle_time_80_u_0.01, mat_bottle_time_80_u_0.01$Generation,
mat_bottle_time_80_u_0.01$delta.var.adm.props,
mat_bottle_time_80_u_0.01$alpha, mat_bottle_time_80_u_0.01$alpha,ligne = T,
"Delta var adm props en fonction des générations\n bottleneck = 80, U = 0.01",
legd = TRUE)
# vec_test_tot = c()
# level = "0.1"
# for (level in levels(mat_bottle_time_80_u_0.01$alpha)) {
# mat_tmp = mat_bottle_time_80_u_0.01[which(mat_bottle_time_80_u_0.01$alpha == level),]
# vec_test = c()
# for (gen in 1:102) {
# vec_test[gen] = (1 - mat_tmp$mean.adm.props[gen])
# }
# if (level == levels(mat_bottle_time_80_u_0.01$alpha)[1]) {
# vec_test_tot = cbind(vec_test, rep(level, length(vec_test)), seq(0,101,1))
# }else{
# vec_test = cbind(vec_test, rep(level, length(vec_test)), seq(0,101,1))
# vec_test_tot = rbind(vec_test_tot, vec_test)
# }
# }
# vec_test_tot[, 1] = as.double(vec_test_tot[, 1])
# colnames(vec_test_tot) = c("stat", "simu", "gen")
# vec_test_tot = as.data.frame(vec_test_tot)
# p = ggplot(vec_test_tot, aes(x = gen, y = stat,group = simu))
#
# p = p + geom_line()
# p = p + theme_classic()
# print(p)
print(difftime(Sys.time(), time.ini))
## Time difference of 5.017781 mins
setwd("../new_methis_bottleneck_50000_snp/")
mat_bottle_10.000 = data.frame.bottleneck(lst_ne0 = c(10000), lst_nef = c(10000),
lst_alpha = c(0.01, 0.1, 0.5, 0.9),
lst_nu = c(0.01, 0.25, 0.49),
lst_bottle = seq(10,90,10))
mat_bottle_time_50_10.000 = mat_bottle_10.000[which(mat_bottle_10.000$time_botl == 50),]
mat_bottle_time_20_10.000 = mat_bottle_10.000[which(mat_bottle_10.000$time_botl == 20),]
mat_bottle_time_80_10.000 = mat_bottle_10.000[which(mat_bottle_10.000$time_botl == 80),]
Moyenne de l’hétérozygotie
boucle_plot_bottle(mat_bottle_10.000, mat_bottle_10.000$mean.het.adm, "Heterozygotie", c(20,50,80),10000)
mini_mat_bottle_10.000 = rbind(mat_bottle_time_20_10.000, mat_bottle_time_50_10.000, mat_bottle_time_80_10.000)
p = plot_stat_gen(mini_mat_bottle_10.000, mini_mat_bottle_10.000$Generation,
mini_mat_bottle_10.000$mean.het.adm, c(),
mini_mat_bottle_10.000$alpha_u, line_t = mini_mat_bottle_10.000$time_botl,
"Moyenne He en fonction des générations\n selon différents bottleneck\n (Ne0 = Nef = 10000)",TRUE,
legd = TRUE)
p = improve_plot_bottle(p =p, vec_abline = c(21, 51, 81), lab_col = "U/alpha",
lab_line = "tpsbott", vec_linetype = c("dotted", "longdash", "solid"))
print(p)
print(difftime(Sys.time(), time.ini))
## Time difference of 5.133401 mins
setwd("../new_methis_bottleneck_50000_snp")
options( "scipen"=100)
#/alpha0.1/Nu0.01/bottle10/Ne100000-XXX/Ne100000-100000/simu_1/
mat_bottle_100.000 = data.frame.bottleneck(lst_ne0 = c(100000), lst_nef = c(100000),
lst_alpha = c(0.001, 0.01, 0.1, 0.5, 0.9),
lst_nu = c(0.01, 0.25, 0.49),
lst_bottle = seq(10,90,10))
#seq(10, 90, 10)
# mat_bottle_time_50_100.000 = mat_bottle_100.000[which(mat_bottle_100.000$time_botl == 50),]
# mat_bottle_time_20_100.000 = mat_bottle_100.000[which(mat_bottle_100.000$time_botl == 20),]
# mat_bottle_time_80_100.000 = mat_bottle_100.000[which(mat_bottle_100.000$time_botl == 80),]
options( "scipen"=0)
boucle_plot_bottle(mat_bottle_100.000, mat_bottle_100.000$mean.het.adm, "Heterozygotie", c(20,50,80),"100 000")
print(difftime(Sys.time(), time.ini))
## Time difference of 5.243455 mins
setwd("../multi_new_methis_pop_size_diff/")
mat_cstt_multi = data.frame.stat(seq_ne = c(100,1000,10000), max_gen = 101, max_simu = 100)
mat_cstt_multi = delta.adm.props(mat_cstt_multi)
mat_cstt_multi_mean = data.frame.mean(mat_cstt_multi, max_gen = 101,
seq_ne = c(100,1000,10000),
col_cstt = c(1:3), col_mean = c(4:66))
mat_cstt_multi_100 = mat_cstt_multi[which(mat_cstt_multi$Generation == 100),]
write_cstt_plot("delta.mean.adm.props", mat_cstt_multi_mean,
"delta_mean_adm/delta_mean_ne_cstt_multi", -0.1, 0.01)
## Saving 7 x 5 in image
write_cstt_plot("delta.var.adm.props", mat_cstt_multi_mean,
"delta_var_adm/delta_var_ne_cstt_multi", -0.2, 0.05)
## Saving 7 x 5 in image
p = plot_stat_gen(mat_cstt_multi_mean, mat_cstt_multi_mean$Generation,
mat_cstt_multi_mean$delta.mean.adm.props, mat_cstt_multi_mean$Ne,
mat_cstt_multi_mean$Ne,ligne = T,
"Delta mean adm props en fonction de s1.0\n selon différentes Ne initiales",
legd = TRUE)
print(p)
p = plot_stat_gen(mat_cstt_multi_100, mat_cstt_multi_100$s1.0,
mat_cstt_multi_100$delta.mean.adm.props, mat_cstt_multi_100$Ne,
mat_cstt_multi_100$Ne,ligne = T,
"Delta mean adm props en fonction de s1.0\n selon différentes Ne initiales",
legd = TRUE, size_point = 2)
p = p +geom_line(linetype = "dotted")
p = p + xlab("s1.0") + ylab("delta mean adm props")
print(p)
setwd("../multi_new_methis_pop_increase_Nu_var/")
seq_nu_multi = sprintf("%.2f", c(0.01, 0.25, 0.49))
seq_ne_ini_multi = c(100)
seq_ne_fin_multi = c(200,400,1000,10000)
seq_combi_multi = as.data.frame(outer(seq_ne_ini_multi, seq_ne_fin_multi, FUN="paste", sep="-"))
seq_combi_multi = as.data.frame.vector(unlist(seq_combi_multi))[[1]]
mat_inc_multi = data.frame.increase.nu(seq_ne_ini = seq_ne_ini_multi, seq_combi = seq_combi_multi,
seq_nu = seq_nu_multi, max_simu = 100)
mat_inc_multi = delta.adm.props(mat_inc_multi)
mat_inc_multi_mean = data.frame.mean.u(df_stat = mat_inc_multi,
max_gen = 101,
seq_combi = seq_combi_multi,
seq_u = c(0.01, 0.25, 0.49),
col_mean = c(4:67,69:70),
col_cstt = c(1:3,68))
mat_inc_multi_100 = mat_inc_multi[which(mat_inc_multi$Generation == 100),]
write_increase_plot("delta.mean.adm.props", mat_inc_multi_mean, "delta_mean_adm/delta_mean_pop_inc_nu_va_multi", -0.1, 0.01)
## Saving 7 x 5 in image
write_increase_plot("delta.var.adm.props", mat_inc_multi_mean, "delta_var_adm/delta_var_pop_inc_nu_va_multi", -0.2, 0.05)
## Saving 7 x 5 in image
p = plot_stat_gen(df = mat_inc_multi_mean, gen = mat_inc_multi_mean$Generation,
stat = mat_inc_multi_mean$delta.mean.adm.props,
color_col = mat_inc_multi_mean$Ne,
line_t = mat_inc_multi_mean$U,ligne = T,
titre = "Delta mean adm props en fonction de s1.0\n selon différentes Ne initiales",
legd = TRUE)
print(p)
p = plot_stat_gen(df = mat_inc_multi_100, gen = mat_inc_multi_100$s1.0,
stat = mat_inc_multi_100$delta.mean.adm.props,
color_col = mat_inc_multi_100$Ne,
line_t = mat_inc_multi_100$U,ligne = T,
titre = "Delta mean adm props en fonction de s1.0\n selon différentes Ne initiales",
legd = TRUE, size_point = 2)
p = p +geom_line(linetype = "dotted")
p = p + xlab("s1.0") + ylab("delta mean adm props")
print(p)
setwd("../multi_new_methis_bottleneck/")
mat_bot_multi = data.frame.bottleneck(lst_ne0 = c(1000), lst_nef = c(1000),
lst_alpha = c(0.1, 0.5),
lst_nu = c(0.01, 0.25, 0.49),
lst_bottle = c(80), max_simu = 100)
mat_bot_multi = delta.adm.props(mat_bot_multi)
mat_bot_multi_mean = data.frame.mean.bottle(df_stat = mat_bot_multi,
max_gen = 101,
seq_combi = c("1000-1000"),
seq_u = c(0.01, 0.25, 0.49),
seq_alpha = c(0.1, 0.5),
seq_bottle = c(80),
col_mean = c(4:67,73:74),
col_cstt = c(1:3, 68:72))
mat_bot_multi_100 = mat_bot_multi[which(mat_bot_multi$Generation == 100),]
write_bottle_plot("delta.mean.adm.props", mat_bot_multi_mean,
"delta_mean_adm/delta_mean_bottleneck_multi", -0.1, 0.01)
## Saving 7 x 5 in image
write_bottle_plot("delta.var.adm.props", mat_bot_multi_mean,
"delta_var_adm/delta_var_bottleneck_multi", -0.2, 0.05)
## Saving 7 x 5 in image
p = plot_stat_gen(df = mat_bot_multi_mean, gen = mat_bot_multi_mean$Generation,
stat = mat_bot_multi_mean$delta.mean.adm.props,
color_col = mat_bot_multi_mean$U,
line_t = mat_bot_multi_mean$alpha,ligne = T,
titre = "Delta mean adm props en fonction de s1.0\n selon différentes Ne initiales",
legd = TRUE)
print(p)
p = plot_stat_gen(df = mat_bot_multi_100, gen = mat_bot_multi_100$s1.0,
stat = mat_bot_multi_100$delta.mean.adm.props,
color_col = mat_bot_multi_100$U,
line_t = mat_bot_multi_100$alpha,ligne = T,
titre = "Delta mean adm props en fonction de s1.0\n selon différentes Ne initiales",
legd = TRUE, size_point = 2)
p = p +geom_line(linetype = "dotted")
p = p + xlab("s1.0") + ylab("delta mean adm props")
print(p)
print(difftime(Sys.time(), time.ini))
## Time difference of 9.3342 mins